home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / convolve.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  10.9 KB  |  442 lines
  1. /*
  2.  * C code from the article
  3.  * "Fast Convolution with Packed Lookup Tables"
  4.  * by George Wolberg and Henry Massalin,
  5.  * (wolcc@cunyvm.cuny.edu and qua@microunity.com)
  6.  * in "Graphics Gems IV", Academic Press, 1994
  7.  *
  8.  *    Compile: cc convolve.c -o convolve
  9.  *    Execute: convolve in.bw kernel out.bw
  10.  */
  11.  
  12. #include <stdio.h>
  13. #include <stdlib.h>
  14.  
  15. typedef unsigned char    uchar;
  16.  
  17. typedef struct {        /* image data structure     */
  18.     int width;        /* image width    (# cols) */
  19.     int height;        /* image height (# rows) */
  20.     uchar *image;        /* pointer to image data */
  21. } imageS, *imageP;
  22.  
  23. typedef struct {        /* packed lut structure        */
  24.     int lut0[256];        /* stage 0 for    5-pt kernel */
  25.     int lut1[256];        /* stage 1 for 11-pt kernel */
  26.     int lut2[256];        /* stage 2 for 17-pt kernel */
  27.     int bias;        /* accumulated stage biases */
  28.     int stages;        /* # of stages used: 1,2,3  */
  29. } lutS, *lutP;
  30.  
  31. /* definitions */
  32. #define MASK        0x3FF
  33. #define ROUNDD        1
  34. #define PACK(A,B,C)    (((A)<<20) + ((B)<<10) + (C))
  35. #define INT(A)        ((int) ((A)*262144+32768) >> 16)
  36. #define CLAMP(A,L,H)    ((A)<=(L) ? (L) : (A)<=(H) ? (A) : (H))
  37. #define ABS(A)        ((A) >= 0 ? (A) : -(A))
  38.  
  39. /* declarations for convolution functions */
  40. void    convolve();
  41. void    initPackedLuts();
  42. void    fastconv();
  43.  
  44. /* declarations for image utility functions */
  45. imageP    allocImage();
  46. imageP    readImage();
  47. int    saveImage();
  48. void    freeImage();
  49.  
  50. /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  51.  * main:
  52.  *
  53.  * Main function to collect input image and kernel values.
  54.  * Pass them to convolve() and save result in output file.
  55.  */
  56. main(argc, argv)
  57. int      argc;
  58. char    **argv;
  59. {
  60.     int    n;
  61.     imageP    I1, I2;
  62.     float    kernel[9];
  63.     char    buf[80];
  64.     FILE    *fp;
  65.  
  66.     /* make sure the user invokes this program properly */
  67.     if(argc != 4) {
  68.         fprintf(stderr, "Usage: convolve in.bw kernel out.bw\n");
  69.         exit(1);
  70.     }
  71.  
  72.     /* read input image */
  73.     if((I1=readImage(argv[1])) == NULL) {
  74.         fprintf(stderr, "Can't read input file %s\n", argv[1]);
  75.         exit(1);
  76.     }
  77.  
  78.     /* read kernel: n lines in file specify a (2n-1)-point kernel
  79.      * Don't exceed 9 kernel values (17-point symmetric kernel is limit)
  80.      */
  81.     if((fp=fopen(argv[2], "r")) == NULL) {
  82.         fprintf(stderr, "Can't read kernel file %s\n", argv[2]);
  83.         exit(1);
  84.     }
  85.     for(n=0; n<9 && fgets(buf, 80, fp); n++) kernel[n] = atof(buf);
  86.  
  87.     /* convolve input I1 with fast convolver */
  88.     I2 = allocImage(I1->width, I1->height);
  89.     convolve(I1, kernel, n, I2);
  90.  
  91.     /* save output to a file */
  92.     if(saveImage(I2, argv[3]) == NULL) {
  93.         fprintf(stderr, "Can't save output file %s\n", argv[3]);
  94.         exit(1);
  95.     }
  96. }
  97.  
  98.  
  99.  
  100. /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  101.  * convolve:
  102.  *
  103.  * Convolve input image I1 with kernel, a (2n-1)-point symmetric filter
  104.  * kernel containing n entries: h[i] = kernel[ |i| ] for -n < i < n.
  105.  * Output is stored in I2.
  106.  */
  107. void
  108. convolve(I1, kernel, n, I2)
  109. imageP     I1, I2;
  110. float    *kernel;
  111. int     n;
  112. {
  113.     int    x, y, w, h;
  114.     uchar    *src, *dst;
  115.     imageP    II;
  116.     lutS    luts;
  117.  
  118.     w = I1->width;                /* image width        */
  119.     h = I1->height;                /* image height        */
  120.  
  121.     II = allocImage(w, h);            /* reserve tmp image    */
  122.     initPackedLuts(kernel, n, &luts);    /* init packed luts    */
  123.  
  124.     for(y=0; y<h; y++) {            /* process all rows    */
  125.         src = I1->image + y*w;        /* ptr to input     row    */
  126.         dst = II->image + y*w;        /* ptr to output row    */
  127.         fastconv(src, w, 1, &luts, dst);/* w pixels; stride=1    */
  128.     }
  129.  
  130.     for(x=0; x<w; x++) {            /* process all columns    */
  131.         src = II->image + x;        /* ptr to input     column */
  132.         dst = I2->image + x;        /* ptr to output column */
  133.         fastconv(src, h, w, &luts, dst);/* h pixels; stride=w    */
  134.     }
  135.  
  136.     freeImage(II);                /* free temporary image */
  137. }
  138.  
  139.  
  140.  
  141. /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  142.  * initPackedLuts:
  143.  *
  144.  * Initialize scaled and packed lookup tables in lut.
  145.  * Permit up to 3 cascaded stages for the following kernel sizes:
  146.  *    stage 0:  5-point kernel
  147.  *    stage 1: 11-point kernel
  148.  *    stage 2: 17-point kernel
  149.  * lut->lut0 <== packed entries (i*k2, i*k1, .5*i*k0), for i in [0, 255]
  150.  * lut->lut1 <== packed entries (i*k5, i*k4,    i*k3), for i in [0, 255]
  151.  * lut->lut2 <== packed entries (i*k8, i*k7,    i*k6), for i in [0, 255]
  152.  * where k0,...k8 are taken in sequence from kernel[].
  153.  *
  154.  * Note that in lut0, k0 is halved since it corresponds to the center
  155.  * pixel's kernel value and it appears in both fwd0 and rev0 (see gem).
  156.  */
  157. static void
  158. initPackedLuts(kernel, n, luts)
  159. float    *kernel;
  160. int     n;
  161. lutP     luts;
  162. {
  163.     int    i, k, s, *lut;
  164.     int    b1, b2, b3;
  165.     float    k1, k2, k3;
  166.     float    sum;
  167.  
  168.     /* enforce flat-field response constraint: sum of kernel values = 1 */
  169.     sum = kernel[0];
  170.     for(i=1; i<n; i++) sum += 2*kernel[i];    /* account for symmetry */
  171.     if(ABS(sum - 1) > .001)
  172.         fprintf(stderr, "Warning: filter sum != 1 (=%f)\n", sum);
  173.  
  174.     /* init bias added to fields to avoid negative numbers (underflow) */
  175.     luts->bias = 0;
  176.  
  177.     /* set up lut stages, 3 kernel values at a time */
  178.     for(k=s=0; k<n; s++) {            /* init lut (stage s)    */
  179.         k1 = (k < n) ? kernel[k++] : 0;
  180.         k2 = (k < n) ? kernel[k++] : 0;
  181.         k3 = (k < n) ? kernel[k++] : 0;
  182.         if(k <= 3) k1 *= .5;        /* kernel[0]: halve k0    */
  183.  
  184.         /* select proper array in lut structure based on stage s */
  185.         switch(s) {
  186.         case 0: lut = luts->lut0;    break;
  187.         case 1: lut = luts->lut1;    break;
  188.         case 2: lut = luts->lut2;    break;
  189.         }
  190.  
  191.         /* check k1,k2,k3 to avoid overflow in 10-bit fields */
  192.         if(ABS(k1) + ABS(k2) + ABS(k3) > 1) {
  193.             fprintf(stderr, "|%f|+|%f|+|%f| > 1\n", k1, k2, k3);
  194.             exit(1);
  195.         }
  196.  
  197.         /* compute bias for each field to avoid underflow */
  198.         b1 = b2 = b3 = 0;
  199.         if(k1 < 0) b1 = -k1 * 1024;
  200.         if(k2 < 0) b2 = -k2 * 1024;
  201.         if(k3 < 0) b3 = -k3 * 1024;
  202.  
  203.         /* luts->bias will be subtracted in convolve() after adding
  204.          * stages; multiply by 2 because of combined effect of fwd
  205.          * and rev terms
  206.          */
  207.         luts->bias += 2*(b1 + b2 + b3);
  208.  
  209.         /* scale and pack kernel values in lut */
  210.         for(i=0; i<256; i++) {
  211.             /*
  212.              * INT(A) forms fixed point field:
  213.              * (A*(1<<18)+(1<<15)) >> 16
  214.              */
  215.             lut[i] = PACK(    INT(i*k3) + b3,
  216.                     INT(i*k2) + b2 + ROUNDD,
  217.                     INT(i*k1) + b1 );
  218.         }
  219.     }
  220.     luts->stages = s;
  221. }
  222.  
  223.  
  224.  
  225. /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  226.  * fastconv:
  227.  *
  228.  * Fast 1D convolver.
  229.  * Convolve len input samples in src with a symmetric kernel packed in luts,
  230.  * a lookup table that is created by initPackedLuts() from kernel values.
  231.  * The output goes into dst.
  232.  */
  233. static void
  234. fastconv(src, len, offst, luts, dst)
  235. int    len, offst;
  236. uchar    *src, *dst;
  237. lutP    luts;
  238. {
  239.     int     x, padlen, val, bias;
  240.     int     fwd0, fwd1, fwd2;
  241.     int     rev0, rev1, rev2;
  242.     int    *lut0, *lut1, *lut2;
  243.     uchar    *p1, *p2, *ip, *op;
  244.     uchar     buf[1024];
  245.  
  246.     /* copy and pad src into buf with padlen elements on each end */
  247.     padlen = 3*(luts->stages) - 1;
  248.     p1 = src;        /* pointer to row (or column) of input    */
  249.     p2 = buf;        /* pointer to row of padded buffer    */
  250.     for(x=0; x<padlen; x++) /* pad left side: replicate first pixel */
  251.         *p2++ = *p1;
  252.     for(x=0; x<len; x++) {    /* copy input row (or column)        */
  253.         *p2++ = *p1;
  254.          p1  +=     offst;
  255.     }
  256.     p1 -= offst;        /* point to last valid input pixel    */
  257.     for(x=0; x<padlen; x++) /* pad right side: replicate last pixel */
  258.         *p2++ = *p1;
  259.  
  260.     /* initialize input and output pointers, ip and op, respectively */
  261.     ip = buf;
  262.     op = dst;
  263.  
  264.     /* bias was added to lut entries to deal with negative kernel values */
  265.     bias = luts->bias;
  266.  
  267.     switch(luts->stages) {
  268.     case 1:        /* 5-pt kernel */
  269.         lut0 = luts->lut0;
  270.  
  271.         ip  += 2;    /* ip[0] is center pixel */
  272.         fwd0 = (lut0[ip[-2]] >> 10) + lut0[ip[-1]];
  273.         rev0 = (lut0[ip[ 0]] << 10) + lut0[ip[ 1]];
  274.  
  275.         while(len--) {
  276.             fwd0 = (fwd0 >> 10) + lut0[ip[0]];
  277.             rev0 = (rev0 << 10) + lut0[ip[2]];
  278.             val = ((fwd0 & MASK) + ((rev0 >> 20) & MASK) - bias)
  279.                 >> 2;
  280.             *op = CLAMP(val, 0, 255);
  281.  
  282.             ip++;
  283.             op += offst;
  284.         }
  285.         break;
  286.     case 2:        /* 11-pt kernel */
  287.         lut0 = luts->lut0;
  288.         lut1 = luts->lut1;
  289.  
  290.         ip  += 5;    /* ip[0] is center pixel */
  291.         fwd0 = (lut0[ip[-2]] >> 10) + lut0[ip[-1]];
  292.         rev0 = (lut0[ip[ 0]] << 10) + lut0[ip[ 1]];
  293.  
  294.         fwd1 = (lut1[ip[-5]] >> 10) + lut1[ip[-4]];
  295.         rev1 = (lut1[ip[ 3]] << 10) + lut1[ip[ 4]];
  296.  
  297.         while(len--) {
  298.             fwd0 = (fwd0 >> 10) + lut0[ip[0]];
  299.             rev0 = (rev0 << 10) + lut0[ip[2]];
  300.  
  301.             fwd1 = (fwd1 >> 10) + lut1[ip[-3]];
  302.             rev1 = (rev1 << 10) + lut1[ip[ 5]];
  303.  
  304.             val  =    ((fwd0 & MASK) + ((rev0 >> 20) & MASK)
  305.                 +(fwd1 & MASK) + ((rev1 >> 20) & MASK) - bias)
  306.                 >> 2;
  307.             *op = CLAMP(val, 0, 255);
  308.  
  309.             ip++;
  310.             op += offst;
  311.         }
  312.         break;
  313.     case 3:        /* 17-pt kernel */
  314.         lut0 = luts->lut0;
  315.         lut1 = luts->lut1;
  316.         lut2 = luts->lut2;
  317.  
  318.         ip  += 8;    /* ip[0] is center pixel */
  319.         fwd0 = (lut0[ip[-2]] >> 10) + lut0[ip[-1]];
  320.         rev0 = (lut0[ip[ 0]] << 10) + lut0[ip[ 1]];
  321.  
  322.         fwd1 = (lut1[ip[-5]] >> 10) + lut1[ip[-4]];
  323.         rev1 = (lut1[ip[ 3]] << 10) + lut1[ip[ 4]];
  324.  
  325.         fwd2 = (lut2[ip[-8]] >> 10) + lut2[ip[-7]];
  326.         rev2 = (lut2[ip[ 6]] << 10) + lut2[ip[ 7]];
  327.  
  328.         while(len--) {
  329.             fwd0 = (fwd0 >> 10) + lut0[ip[0]];
  330.             rev0 = (rev0 << 10) + lut0[ip[2]];
  331.  
  332.             fwd1 = (fwd1 >> 10) + lut1[ip[-3]];
  333.             rev1 = (rev1 << 10) + lut1[ip[ 5]];
  334.  
  335.             fwd2 = (fwd2 >> 10) + lut2[ip[-6]];
  336.             rev2 = (rev2 << 10) + lut2[ip[ 8]];
  337.  
  338.             val  =    ((fwd0 & MASK) + ((rev0 >> 20) & MASK)
  339.                 +(fwd1 & MASK) + ((rev1 >> 20) & MASK)
  340.                 +(fwd2 & MASK) + ((rev2 >> 20) & MASK) - bias)
  341.                 >> 2;
  342.             *op = CLAMP(val, 0, 255);
  343.  
  344.             ip++;
  345.             op += offst;
  346.         }
  347.         break;
  348.     }
  349. }
  350.  
  351.  
  352.  
  353. /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  354.  * readImage:
  355.  *
  356.  * Read an image from file.
  357.  * Format: two integers to specify width and height, followed by uchar data.
  358.  * Return image structure pointer.
  359.  */
  360. imageP
  361. readImage(file)
  362. char    *file;
  363. {
  364.     int     sz[2];
  365.     FILE    *fp;
  366.     imageP     I = NULL;
  367.  
  368.     /* open file for reading */
  369.     if((fp = fopen(file, "r")) != NULL) {    /* open file for read    */
  370.         fread(sz, sizeof(int), 2, fp);    /* read image dimensions*/
  371.         I = allocImage( sz[0],sz[1]);    /* init image structure */
  372.         fread(I->image, sz[0],sz[1],fp);/* read data into I    */
  373.         fclose(fp);            /* close image file    */
  374.     }
  375.     return(I);                /* return image pointer */
  376. }
  377.  
  378.  
  379.  
  380. /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  381.  * saveImage:
  382.  *
  383.  * Save image I into file.
  384.  * Return NULL for failure, 1 for success.
  385.  */
  386. int
  387. saveImage(I, file)
  388. imageP     I;
  389. char    *file;
  390. {
  391.     int     sz[2], status = NULL;
  392.     FILE    *fp;
  393.  
  394.     if((fp = fopen(file, "w")) != NULL) {    /* open file for save    */
  395.         sz[0] = I->width;
  396.         sz[1] = I->height;
  397.         fwrite(sz, sizeof(int), 2, fp); /* write dimensions    */
  398.         fwrite(I->image,sz[0],sz[1],fp);/* write image data    */
  399.         fclose(fp);            /* close image file    */
  400.         status = 1;            /* register success    */
  401.     }
  402.     return(status);
  403. }
  404.  
  405.  
  406.  
  407. /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  408.  * allocImage:
  409.  *
  410.  * Allocate space for an uchar image of width w and height h.
  411.  * Return image structure pointer.
  412.  */
  413. imageP
  414. allocImage(w, h)
  415. int w, h;
  416. {
  417.     imageP     I;
  418.  
  419.     /* allocate memory for image data structure */
  420.     if((I = (imageP) malloc(sizeof(imageS))) != NULL) {
  421.         I->width  = w;            /* init width        */
  422.         I->height = h;            /* init height        */
  423.         I->image  =(uchar*) malloc(w*h);/* init data pointer    */
  424.     }
  425.     return(I);                /* return image pointer */
  426. }
  427.  
  428.  
  429.  
  430. /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  431.  * freeImage:
  432.  *
  433.  * Free image memory.
  434.  */
  435. void
  436. freeImage(I)
  437. imageP I;
  438. {
  439.     free((char *) I->image);
  440.     free((char *) I);
  441. }
  442.